function [Ypred,numsv] = dqsvm(X, Y, ktype, C)
% SVM function [Ypred] = dqsvm(data, label, ktype, C)
%       classifies the unlabeled instances in data according
%       to the labeled ones using support vector machine
%
%       label contains the class for first part of data
%       and Ypred contains the results for the rest
%
%       ktype == 1 linear
%                     2 polynomial
%                     3 rbf
%
%       version 0.1, 04/06/2014
%       Suqi Liu

normX = sqrt(sum(X.^2,2));
X = bsxfun(@rdivide,X,normX);

l = size(Y,1);
n = size(X,1);

Xtr = X(1:l,:);
Xts = X(l+1:n,:);

dim = 3;
sigma = 1;

if (ktype ==1)
    K = Xtr*Xtr';
end
if (ktype == 2)
    K = (1+Xtr*Xtr').^dim;
end
if (ktype == 3)
    Ktmp = Xtr*Xtr';
    xs = diag(Ktmp);
    K = exp((2*Ktmp - repmat(xs,1,l) - repmat(xs',l,1))/(2*sigma));
    clear Ktmp;
end

Q = (Y*Y').*K;
clear K;

TolVal = 1e-8;

tic
if (nargin == 3)
    a = quadprog(Q,-ones(l,1),[],[],Y',0,zeros(l,1),[]);
    svi = find(a>100*sqrt(TolVal));
else
    a = quadprog(Q,-ones(l,1),[],[],Y',0,zeros(l,1),C*ones(l,1));
    svi = find(a<C-100*sqrt(TolVal) & a>100*sqrt(TolVal));
end
optime = toc

numsv = length(svi);

b = (1-a'*Q(:,svi))*Y(svi)/numsv;

clear Q;

if (ktype ==1)
    T = Xts*Xtr';
end
if (ktype == 2)
    T = (1+Xts*Xtr').^dim;
end
if (ktype == 3)
    Ttmp = Xts*Xtr';
    ts = sum(Xts.*Xts,2);
    T = exp((2*Ttmp - repmat(ts,1,l) - repmat(xs',n-l,1))/(2*sigma));
    clear Ttmp;
    clear ts;
end

Ypred = sign(T*(Y.*a)+repmat(b,n-l,1));
end
